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A FORTRAN  PROGRAM  TO  FIT  A COMPLEX- VALUED  TRANSFER  FUNCTION 


INTRODUCTION 

A method  Is  programmed  to  fit  a transfer  function  with  up  to  eight 
parameters:  gain,  delay,  3 leads  (or  zeros),  and  3 lags  (or  poles).  A 
modified  version  of  a Fortran  program  by  Donald  W.  Marquardt  Cl)  for 
least-squares  estimation  of  nonlinear  parameters  Is  the  basic  program 
used.  The  program  requires  a subroutine  to  define  the  residual  and  the 
partial  derivatives  of  the  function  to  be  fitted  and  to  evaluate  these 
real- valued  quantities  at  each  of  N points.  This  paper  discusses  the 
assumptions  used  and  the  structure  of  the  subroutine  to  fit  a complex 
valued  function  with  a program  written  to  fit  real  valued  functions. 


ASSUMPTIONS  MADE  IN  THE  METHOD 

The  norm  to  be  minimized  is: 

$ = JjlH(jwp  - H(j“i)l^  (1) 

where  w Is  the  frequency  in  radians,  H is  the  observed  value  of  the 
transfer  function,  and  ft  is  the  predicted  value  of  the  transfer  function 
using  the  estimated  values  of  the  parameters. 

The  observed  values  can  be  written  as: 

H(j  0)^)  * Ui  + jVi  (2) 

and  the  predicted  values  can  be  written  as: 

H(ju)i)  = Ri(0)  jQi(e)  (3) 

where  0 represents  the  vector  of  parameters  and  the  R^(0)  and  Q (0)  are 
real  functions  of  these  parameters  and  the  frequency.  These  functions 
are  given  in  Appendix  A for  the  eight-parameter  transfer  function  con- 
sidered herein.  Using  equations  2 and  3,  equation  1 can  be  written  as: 

^ ",  f^l  “ ^ tv.  - Q,(0)]2  (4) 

1=)  1 1 i=l  1 1 

The  main  program  requires  that  the  sum  of  squares  of  residuals  to  be 
minimized  be  of  the  form 


1.  Marquardt,  D.  W.  An  algorithm  for  least-squares  estimation  of 
nonlinear  parameters.  J Soc  Indust  Appl  Math  II  (No.  2) : 431-441  (1963)  . 
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(5) 


Z Res^  = E (observed  - function)^  . 

In  equation  4,  the  function  is  (6)  for  the  real  part  and  (6)  for 
the  imaginary  part.  These  two  functions  require  different  computations 
for  the  residual  and  different  sets  of  partial  derivatives  to  be 
evaluated.  These  functions  and  their  partial  derivatives  with  respect 
to  the  parameters  are  given  In  Appendix  A for  the  eight-parameter  trans- 
fer function. 

Under  the  assumption  that  the  real  and  imaginary  parts  of  H(ja)j^) 
are  Independent  for  each  and  that  the  are  independent  for 

different  u),  the  problem  of  fitting  a complex-valued  function  to  n fre- 
quencies can  be  considered  as  a problem  of  fitting  one  of  two  real- 
valued functions  to  N > 2n  Independent  data  points.  When  only  ensemble 
averaging  Is  done  to  calculate  the  observed  values  of  the  transfer 
function,  the  approximate  Independence  of  the  values  at  the  Fourier  fre- 
quencies Is  an  imnedlate  extension  of  the  approximate  Independence  for 
the  perlodogram.  When  tapering  or  smoothing  is  done,  this  will  no 
f longer  be  true,  except  for  smoothed  estimates  at  frequencies  not  using 

any  common  unsmoothed  value.  Therefore,  ensemble  averaging  Is  the 
[ recommended  method  of  obtaining  the  observed  values.  When  only  a sdngle 

sample  Is  available,  the  use  of  the  original  Bartlett.. wlnd.pw  that  Ip- 
r volves  dividing  the  sample  Into  several  subsamples  of  equa^  length  and 

averaging  the  perlodogram  estimates  from  these  subsamples  Is  tb^e 
recommended  alternative.  

The  use  of  averaging  still  does  not  assure  any  reasonability  of 
I the  assumption  of  Independence  of  the  real  and  Imaginary  parts  at  the 

I same  frequency.  This  appears  to  be  analytically  Intractable  to’show, 

I but  preliminary  results  from  some  Monte  Carlo  simulation  Indicate  that 

|.  the  assumption  Is  probably  not  unreasonable. 

Finally,  the  reasonableness  of  using  equal  weights,  at  all*  fre- 
f quencles  should  be  considered.  Assuming  that  the  variances  of  the  real 

I and  Imaginary  parts  of  all  frequencies  are  equal,  then  equal  weights 

should  be  used.  Demonstrating  equality  of  variances  appears  to  be 
analytically  Intractable,  although  preliminary  results  from  the  simula- 
tion data  suggest  the  assumption  may  be  reasonably  sound.  Therefore, 
estimating  the  parameters  of  a transfer  function  by  minimizing  equation 
4 and  considering  the  2n  > N data  values  as  Independent  with  equal 
variances  seems  to  be  reasonable. 


STRUCTURE  OF  THE  SUBROUTINE 

The  required  standard  Information  available  to  the  subroutine  con- 
sists of  the  vector  of  parameter  estimates  from  the  previous  iteration, 
the  vector  of  N values  of  the  dependent  variable  (Y) , the  matrix  of  N 
I values  of  one  or  more  Independent  variables  (X) , and  the  Index  of  the 


2 


r 


point;  whose  residual  and  partial  derivatives  are  currently  being 
evaluated.  The  additional  infomatlon  needed  for  this  estimation 
problem  Is  whether  the  current  value  Is  a real  part  or  imaginary  part 
of  an  observed  value.  This  Is  obtained  by  coding  an  extra  Independent 
variable,  1 If  real  and  2 if  imaginary.  Thus  each  of  the  N data 
vectors  is  made  up  of  three  quantities:  the  real  or  Imaginary  part 
of  an  observed  value,  the  frequency  In  radians  associated  with  the  ob- 
served value,  and  an  Indicator  of  whether  the  data  vector  Is  for  a real 
part  or  an  imaginary  part  (coded  1 or  2) . 


Since  the  real  and  imaginary  parts  of  the  residuals  and  the  partial 
derivatives  are  determined  separately,  the  computations  in  the  sub- 
routine are  done  In  two  sections,  one  for  the  real  part  and  the  other 
for  the  imaginary  part.  Therefore,  the  subroutine  Is  a two-function 
subroutine  with  an  Indicator  to  specify  which  function  to  use.  An 
example  of  a subroutine  which  incorporates  the  Ideas  described  above 
for  the  eight-parameter  transfer  function  Is  shown  In  Appendix  B. 
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APPENDIX  A 


FUNCTIONS  OF  THE  REAL  AND  IMAGINARY  PARTS  OF  THE 
EIGHT-PARAMETER  TRANSFER  FUNCTION  AND  THEIR  PARTIAL 
DERIVATIVES  WITH  RESPECT  TO  THE  PARAMETERS 

The  transfer  function  can  be  written  as: 

M(Ja.i)  - K[H-Y,au>^)  -I-  Yg  (ju)^)^] 


[l  + o^(Ja)i;  + 


Writing  the  exponential  as: 

g jTtOl  ^ _ J ^ 

multiplying  the  numerator  and  denominator  by  the  conjugate  of  the  de- 
nominator, and  collecting  terms  gives: 

H(ja)i)»  Ri  + jQi 

where,  dropping  the  1 subscript, 

R - K[(AB  + CD)  cos  TO)  + (AD-BC)  sin  TuJ/Z 
Q = K[(AD  - BC)  cos  TO)  - (AB  + CD)  sin  tojJ/Z 


A • 1 ~ 

B * 1 - 

C ■ (oCoj  “ 0301^) 

D - oj(yi  - Y3“^) 

2 2 
Z - A + C 

Denoting  the  partial  derivatives  of  R and  Q with  respect  to  the 
parameter  9 as  r|0  and  q|6  and  using  the  above  defined  terms,  the 
partial  derivatives  are: 
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rIk  = R/K 

Rh  = Ku)t(AD-  BC)  cos  TO)  - (AB  + CD)  sin  tidJ/Z 

rIy^*  KwtC  cos  TU)  + A sin  Ta)]/Z 

rIy^*  - A cos  TO)  + C sin  t(u]/Z 

rIy  = Kt»)^  [ - C cos  TU)  - A sin  tul/Z 
3 

R(a^  =.(0  [K(D  cos  TU)  - B sin  tu))  - 2CR]/Z 

r|o^=  u)^  [K  ( - B cos  TU)  - D sin  tw)  + 2AR]/Z 

R|a  = u)^  [K(  - D cos  tu)  + B sin  tu))  + 2CR]/Z 
3 

and 

QlK  = Q/K 

qIt  = Ku)t  - (AB  + CD)  cos  tu)  - (AD  - BC)  sin  TU)]/Z 

Q|y  = K'i)[A  cos  XU)  - C sin  tu)]/Z 

q(y  = Ku)2[  C cos  TU)  + A sin  TU)]/Z 
2 

qIy  * - A cos  TU)  + C sin  tu)]/z 

3 

io[K(-B  cos  tu)  - D sin  Tio)  - 2CQ]/Z 

qIo  = u)2 [K(-D  cos  TU)  + B sin  too)  + 2AQ]/Z 
2 

Q|a  ■ u)3[K(B  cos  tu)  + D sin  too)  + 2CQ]/Z 
3 
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APPENDIX  B 


LISTING  OF  SUBROUTINE 

The  following  subroutine  Is  used  In  conjunction  with  our  modified 
version  of  a program  for  least-squares  estimation  of  nonlinear  param- 
eters. The  section  for  the  first  key  puts  out  the  labeling  for  the 
full  model  form  and  sets  up  the  test  value  for  the  Indicator  of  real  or 
Imaginary  part.  The  second  key  section  estimates  the  function  and  the 
residual  and  the  third  key  section  estimates  the  partial  derivatives. 

SUBROUTINE  MODEL (KEY) 

C FITTING  COMPLEX  VALUED  TRANSFER  FUNCTION:  EQUAL  WEIGHTS: 

IMPLICIT  REAL*8(A-H,0-Z) 

COMMON  X(300 ,10) ,Y(300) ,B(30) ,P(30) ,PRNT(5) ,HDR2(18) ,F,RES,NI(40) , 
1I,N,K,IP,NPRNT 
DIMENSION  XHDR(18) 

DATA  XHDR/14AHM0DEL  FOR  FITTING  IN  THE  FREQUENCY  DOMAIN ;S“J  OMEGA 
1 H(S)-K(l+€(l)S+...+G(3)S**3)EXP(-TS)  / (1+A(1) 

2S+. . .+A(3)S**3)  / 

GO  TO  (10,20,30),  KEY 
10  DO  5 1=1,18 
5 HDR2(I)=XHDR(I) 

REAL=1.  DO 
RETURN 
C B(1)=K 
C B(2)-T 
C B(3)=G(1) 

C B(4)-G(2) 
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C B(5)=G(3) 

C B(6)*A(1) 

C B(7)-A(2) 

C B(8)*A(3) 

20  W2-X(I,1)**2 
•ni=B(2)*X(I,l) 

A-1.D0-B(7)*W2 

E-1.D0-B(4)*W2 

C=X(1,1)*(B(6)-B(8)*W2) 

D-X(I,1)*(B(3)-B(5)*W2) 

Z-A**2+G**2 

CS-DCOS(TW) 

SN-DSIN(TW) 

ABCD-A*E+C*D 

ADBC-A*D-E*C 

IF(X(I,2).GT.REAL)  GO  TO  25 
F=B (1) * (ABCD*CS+ADBC*SN) /Z 
RES-(Y(I)-F) 

RETURN 

25  F-B(1)*(ADBC*CS-ABCD*SN)/Z 
RES-(Y(I)-F) 

RETURN 
30  W-X(I,1) 

W3-W*W2 

PK-B(l) 


IF(X(I,2).GT.REAL)  GO  TO  35 
P(1)=F/PK 


P (2)=PK*W* (ADBC*CS-ABCD*SN) /Z 
P(3)=PK*W*(C*CS+A*SN)/Z 
P (4 ) “PK*W2* (C*SN-A*CS) /Z 
P (5)=PK*W3* (-C*CS-A*SN) /Z 
P (6)=W* (PK* (D*CS-E*SN)-2 . D0*C*F) /Z 
P ( 7 ) =W2* (PK* (- E*CS-D*SN)+2 . D0*A*F) /Z 
P (8)  =»W3*  (PK*  (-D*CS+E*SN)+2 . DO*C*F)  /Z 
RETURN 

35  P(1)-F/PK 

P (2) *PK*W* (-ADBC*SN-ABCD*CS) /Z 

P(3)-PK*W*(A*CS-C*SN) /Z 

P(4)=PK*W2*(C*CS+A*SN)/Z 

P(5)-PK*W3*(-A*CS+C*SN)/Z 

P (6 ) -W* (PK* (- E*CS-D*SN ) - 2 . DO*C*F) /Z 

P (7)=W2* (PK* (-D*CS+E*SN)+2 . DO*A*F) /Z 

P(8)-W3*(PK*(E*CS  D*SN)+2.D0*C*F)/Z 

RETURN 

END 


